function fval = int_capital(L, a_b, a_s, EPS_L_BOUND, EPS_H_BOUND, S_PARAM, M_PARAM, THETA, CHI_PI, CHI_PI_k)

sig = S_PARAM;
mu  = M_PARAM;
f   = @(x, a) 1 ./ (x * sig * sqrt(2 * 3.14159265359))...
    .* exp(- (log(x) - mu - a).^2 / 2 / sig^2)...
    .* (x >= EPS_L_BOUND) .* (x <= EPS_H_BOUND);
int_method = 'auto';
Rel_Tol    = 1e-5;
Abs_Tol    = 1e-6;


fun.l_bound      = @(a) EPS_L_BOUND; % lower bound
fun.h_bound      = @(a) EPS_H_BOUND; % upper bound
fun.threshold    = @(x, L, CHI_PI, CHI_PI_k) L ./ (1 - THETA - CHI_PI - CHI_PI_k) - THETA ./ (1 - THETA - CHI_PI - CHI_PI_k) .* x; 


fval = integral2(@(eps, eps_t) (eps_t - eps) .* f(eps, a_s) .* f(eps_t, a_b),...
                    fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
                    @(eps)eps, @(eps)fun.threshold(eps, L, CHI_PI, CHI_PI_k),...
                    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);


end